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I. Introduction 

Characteristic-based explicit and implicit total variation diminishing (TVD) schemes for the two- 
dimensional compressible Euler equations have recently been developed [1,2]. This is a generalization of 
recent work of Roe and Davis [3,4] to a wider class of symmetric (non-upwind) TVD schemes other than 
Lax-Wendroff. Roe and Davis’s schemes can be viewed as a subset of the class of explicit methods. The 
main properties of the present class of schemes are that they can be implicit, and, when steady-state cal- 
culations are sought, the numerical solution is independent of the time step. In reference [2], a comparison 
of a linearized form of the present implicit symmetric TVD scheme with an implicit upwind TVD scheme 
originally developed by Harten [5,6] and modified by Yee [7] was given. The results favored the symmetric 
method. It was found that the symmetric method is just as accurate as the upwind method while requiring 
less computational effort. It is emphasized that the generic use of the notion upwind and symmetric TVD 
schemes here pertains to the schemes without the limiter present. With a limiter present, an upwind TVD 
scheme no longer has its traditional upwinding meaning. The same situation also applies to symmetric TVD 
schemes. Another way of distinguishing an upwind from a symmetric TVD scheme is that the numerical 
dissipation term corresponding to an upwind TVD scheme is upwind weighted [5-10], as opposed to the 
numerical dissipation term corresponding to a symmetric TVD scheme which is centered [1-4]. 

Currently, more numerical experiments are being conducted on time-accurate calculations and on the 
effect of grid topology, numerical boundary condition procedures, and different flow conditions on the be- 
havior of the method for steady-state applications. The purpose of this paper is to report our experiences 
with this type of scheme and give some basic guidelines for using the scheme. A description of upwind and 
symmetric TVD schemes, including formulation and extension to the two-dimensional Euler equations of 
gas dynamics in curvilinear coordinates, can be found in references [2,7,11]. This paper contains a brief 
description of the TVD algorithm and a summary of the numerical computations. 


II. Description of Algorithm in One Dimension 

Consider a one-dimensional system of hyperbolic conservation laws 

dU dF(U) 

at dx 


(i) 


Here U and F(U) are column vectors of m components and A — dF/dU . Let the eigenvalues of A be 
(a 1 , a 2 , ...,a m ). Denote R as the matrix whose columns are eigenvectors of A, and R -1 as the inverse of R. 
Let the grid spacing be denoted by Ax such that x — j’Ax, and let Uj+± denotes some symmetric average 

of Uj and {/,■+ 1 (for example, U + 1 = \(U 3 + i + U 3 ) ). Let a 1 . A , R , i, R~\ denote the quantities a\ R, 

R -1 related to A evaluated at i . Define 

Jt 2 

a i+i = R J+i( U i^- U i) ( 2 ) 


as the difference of the local characteristic variables. 

A one-parameter family of conservative explicit and implicit second-order TVD schemes can be written 
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where A ~ , At is the time step and 0 is a parameter. When 0 ^ 0, the scheme is implicit. The numerical 

flux function ± can be expressed as 


J5T i+ i = - 

2 2 


(3b) 


The elements of the denoted by (^- +A ) 5 for a general second-order symmetric TVD scheme [l] are 


t^+i) - -A^(a‘ + i) <?' + i - V»K + i) [«' + 1 - <9' + i 


(4a) 


where a( + , are elements of (2). When 3 =■ 1, the method is best suited for time-accurate calculations, and 
when /? = 0, it is mainly for steady-state applications. The function xp is 


« H z ) 


\z\ \z\ > e 

(z 2 + e 2 )/2e \z\ < c ‘ 


(4b) 


Here tp(z) in (4b) is an entropy correction to \z\ where t is a small positive parameter (see reference [12] for 
a formula for e). Examples of the ‘limiter’ function Q l ± can be expressed as 


Q l j+i = min mod(a'_, ,aj. + i ) + minmod (a‘. + i,a' + § ) - aj. + i 
Q l j+ i =minmod(a'_ i ,a‘ + i) a' + f ) 



= minmod 


2a ‘ - \ - 2a ^+ h > 2a ' + f’iK-i + a ^+ § ) 


(4c) 

(4d) 

(4e) 


The minmod function of a list of arguments is equal to the smallest number in absolute value if the list of 
arguments is of the same sign, or is equal to zero if any arguments are of opposite sign. 

The elements of the 3^ + 1 denoted by (<£*._,_ i )* 7 for a second-order upwind TVD scheme, originally 
developed by Harten, and later modified and generalized by the author [1,13], are 


K+i) 17 = *K+a)(s!+i + *}) - V'K+i + § K+f 

The function a(z) = \${z) + (0 - |)A/?z 2 and 


7L1 = *«+*) 




Examples of the ‘limiter’ function <7* can be expressed as 
y' = minmod(a'_i,a( + i) 

3' = S •max|o,min(2|a' + i|,5 • a'_i),min(|a' + i |,25 • a' _i) |; 5 = sgn(a‘ + i). 


(5a) 

(5b) 

(5c) 

(5d) 

(5e) 


The method of extending scalar TVD schemes to nonlinear systems of hyperbolic conservation laws (1) 
in equation (3) is sometimes referred as the local characteristic approach and is a variant of Roe’s linear 
Riemann solver [14]. The advantages of this approach as opposed to Davis’s simplified approach [4] to 
systems are that: (a) this approach in effect uses scalar schemes on each characteristic field; (b) the limiter 
used need not be the same for each field; e.g., one can use a more compressive limiter for the linear fields 
and use a less compressive limiter for the nonlinear fields; (c) one can even use different schemes for different 
fields; (d) it is more efficient than the exact or approximate Riemann solvers; (e) it provides a natural 
way to linearize the implicit TVD schemes. For the one-dimensional Euler equation of gas dynamics, the 
characteristic fields consist of two nonlinear fields u± c and a linear field u. 



For two-dimensional time-accurate calculations, the explicit scheme f) — 1,0 = 0 is implemented by the 
Strang type of time splitting method [15]. For steady-state numerical study, the implicit scheme /? = 0,0 = 1 
considered here is implemented in a conservative noniterative ADI form [2], The numerical solution is 
independent of the time step. The implicit operator has a regular block tridiagonal structure and the 
resulting block tridiagonal matrix is diagonally dominant. One can modify a classical ADI central difference 
code by simply changing the linear numerical dissipation term into the one designed for the TVD scheme; 
i.e., the third term of equation (3b). See references [7,11] for more details. 

Let the grid indexing at the left boundary in the x-direction be j=0. For the explicit operator, this is a 
five-point scheme in each coordinate direction and one needs values of g l 0 as well as Uo- For simplicity and 
illustration purposes, the values of Uq are assumed to be updated by the procedure described in reference 
[11]. The following three ways of obtaining at the boundaries are discussed. 

0o = °; 9o = 9v 9o = ol 1 i. (6a,b,c) 
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Fig. 1 The Density Distribution (- exact, o com- 
puted) of the shock tube problem using 
limiter (4c-4e) for the symmetric scheme 
and using (5d) for the upwind scheme (top 
to bottom) 


One can interpret (6a) as a first-order TVD 
scheme at the boundary, whereas (6b) is obtained 
by zeroth-order extrapolation. Boundary scheme 
(6c) is the least dissipative among the three. It 
is obtained by simply removing the logic in the 
limiter and using the next availiable a z . 

III. Numerical Results 

3.1 Time Accurate Calculations: The nu- 

merical experiments were mainly performed on 
the one-dimensional shock tube problem and the 
NACA0012 and NACA0018 airfoils. Limiter (4e) 
produces slightly sharper shocks and contact dis- 
continuities than (4c) or (4d). 



X 


Fig. 2 The 299 x 79 C grid for the NACA0018 
airfoil. 


One Dimension : Figure 1 shows the density distribution for the same shock tube problem as used in reference 
[5]. The solid lines are the exact solutions and the circles are the results computed by using limiters (4c-4e) 
for the symmetric TVD scheme and (5d) for the upwind scheme respectively (from top to bottom). With 
the use of the same limiter, the study shows that in general the symmetric scheme is slightly more diffusive 
than the upwind scheme. However, the resolution (not shown) of the symmetric scheme with limiter (4e) 
is better than the upwind scheme with limiter (5c). Note that not all of the limiters that are designed for 
upwind schemes are applicable for symmetric schemes. 

Two Dimensions : In the experiment, a time dependent curved shock wave collides with the NACA0018 
airfoil at an angle of attack a = 30°. The Mach number Ms is approximately 1.5 at the moment of collision. 
The experiment was performed by Mandella and Bershader [ 16] of Stanford University. In the numerical 
simulations, a constant planar incident shock of Ms = 1.5 and a = 30° is used as an initial condition just 
before the leading edge of the airfoil. This work was done jointly with Young Moon of U.C. Berkeley. Figure 
2 shows the 299 x 79 C grid for the NACA0018 airfoil. Figure 3 shows the numerical results in the form of 
density contour plots using (4e) with the symmetric scheme, and (5d) for the nonlinear fields and (5e) for the 
linear fields with the upwind scheme at (approximately) the four sequential instances of the interferograms. 
The incident and reflected shocks, Mach stems, and slip surface on the lower surface are captured within 3 
grid points. Also the vortices at the upper nose and the trailing edge of the airfoil are well captured in the 
simulation. On the upper surface the experimental results show a weak slip surface which it is not captured 
in the simulation. By increasing the grid resolution around that region, this slip surface is also expected to 
be captured. Even though a planar shock was used to model the curved shock, the shape and location of 
the discontinuities compare favorably with the experiments. The shock resolution of the symmetric TVD 
scheme is similar to the upwind scheme except that it is slightly more diffusive. This leads us to believe that 
even closer agreement between the numerical simulation and the experiments is possible if a time dependent 
curved shock is used. A detail description of the physical problem and a more closely simulated result with 
the exact experiment can be found in a paper by Moon and Yee [17]. Aside from the difference in CPU 
time (at least 35 - 40% less), one advantage of symmetric TVD schemes over upwind schemes is that the 
symmetric TVD schemes are less sensitive to numerical boundary condition treatments for higher Mach 
number flows. Figure 4 shows the computations of the symmetric TVD scheme (4e) at a CFL number of 
0.98 for a planar incident shock of Ms = 20 and a = 30°. With the same numerical boundary conditions, the 
upwind scheme diverges; i.e., a characteristic type of boundary condition is needed for the upwind scheme. 

3,2 Steady-State Calculations: Some steady-state computations have been carried out to illustrate 

the sensitivity of the implicit symmetric (and upwind) scheme on the various numerical boundary condition 
procedures, grid topologies and flow conditions. A typical solution computed by the symmetric scheme using 
(4e) is shown in figure 5. Figure 6 shows the 163 x 49 C grid with the outer boundary 24 chord lengths away 
from the body. Some of the numerical study using (4c) and (4d) with (3 = 0 and 0 = 1 using the same grid 
is summarized below. 

Accuracy : To demonstrate how numerical boundary conditions affect the accuracy of the scheme, the compu- 
tations of the NACA0012 airfoil at a freestream Mach number — 0.8 and a — 1.25° with (4c) are shown 
in figure 7. This figure compares the Mach contours of boundary schemes (6a) with (6b). One can see a 
drastic difference in resolution near the solid body. Away from the body, the resolution is indistinguishable. 
A comparison was also made between boundary schemes (6b) and (6c). No visible difference in accuracy 
was found. 

Stability : Computation of the NACA0012 airfoil with freestream Mach numbers ranging from 0.8 to 1.8 with 
and without lift show that boundary schemes (6a) and (6b) have a similar stability and convergence rate. 
However boundary scheme (6c) has a lower stability bound than (6a) and (6b). This is due to the fact that 
boundary scheme (6c) is the least dissipative among the three procedures. 

Grid Topology : The grid shown in figure 1 was generated by a hyperbolic grid method. This grid is fairly 
regular and for the most part is nearly orthogonal. The present method appears to have no difficulty in 
obtaining a fairly accurate solution. If one were to use a highly irregular and nonorthogonal grid, the present 
scheme would require a long time to reach steady state or, in an extreme case, fail to converge. 

Convergence Rate : A numerical study was also conducted on the effect of a freestream Mach number on the 
convergence rate of the scheme. The following is based on 0.7 < M^ < 1.8, 0° < a < 7°, and grid sizes 
ranging from 163 x 33 to 249 x 49 with either a C or O grid. For 0.7 < < 0.8, a Z/ 2 -norm residual of 


i 



ORIGINAL PAGE IS 
OF POOR QUALITY 






7<t 





V "t 


Fig. 3 The Density contours by the upwind scheme (5a,d,e) (middle) and by the 
symmetric scheme (4a, e) (right) for the NACA0018 airfoil with Ms = 1.5, 
a — 30° compared with the interferograms at approximately the same se- 
quential instances. 
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Fig. 4 The density contours (left) and Mach con- 
tours (right) by the symmetric scheme 
(4a, e) for the NACA0018 airfoil with Ms — 
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Fig. 7 The Mach contours computed by boundary schemes (6a) (left) and (6b) 
(right) for the NACA0012 airfoil with M ^ — 0.8, a = 1.25°. 




10 7 can be reached in 600-1200 steps depending on the grid size, time step and a. For 0.85 < < 0.98, 

a residual of 10~ 7 can be reached in 1500-3500 steps. For 1.2 < < 1.8, a residual of 10 7 can be reached 

in 200-500 steps. 


IV. Concluding Remarks 

For problems containing shocks only, numerical experiments for steady-state calculations show that the 
implicit symmetric TVD schemes, while requiring less computational effort, are just as accurate as the 
implicit upwind TVD scheme originally developed by Harten [5] and modified by Yee [7]. However, for 
problems containing both shocks and contact discontinuities, numerical experiments show that the explicit 
symmetric TVD schemes are slightly more diffusive than the explicit upwind TVD scheme especially at the 
contact surfaces. The current study on the various steady-state airfoil calculations also indicates the degree of 
sensitivity of the method to grid topology, freestream Mach number, and angle of attack. Moreover, proper 
numerical boundary condition procedures are essential for numerical simulation of fluid flow. Improper 
treatment of boundary conditions can lead to instability and inaccuracy, even if one starts with a stable 
high-resolution scheme that catered to problems with shocks. Furthermore, for certain flow regimes, different 
linearization or relaxation procedures other than ADI should be investigated to improve the convergence 
rate of steady-state calculations. 


References 

1 H.C. Yee, NASA-TM 86775, July 1985, also to appear, J. Comput. Phys., 1986. 

2 H.C. Yee, Proc. 6th GAMM Conf. Numer. Meth. in Fluid Mechanics, W. Germany, Sept. 1985. 

3 P.L. Roe, ICASE Report No. 84-53, October 1984. 

4 S.F. Davis, ICASE Report No. 84-20, June 1984. 

5 A. Harten, NYU Report, Oct., 1982; SIAM J. Num. Anal, Vol. 21, 1984, pp. 1-23. 

6 A. Harten, J. Comput. Phys., Vol. 49, 1983, pp 357-393. 

7 H.C. Yee, Int. J. Comput. Math. Appl., in press. 

8 P.L. Roe, AMS publications, Lectures in Appl. Math., Vol. 22, 1985. 

9 P.K. Sweby, SIAM J. Num. Analy., Vol. 21, 1984. 

10 S. Osher and S. Chakravarthy, SIAM J. Num. Analy., Vol. 21, 1984. 

11 H.C. Yee and A. Harten, AIAA Paper No. 85-1513, July 1985, to appear AIAA J. 

12 A. Harten and J.M. Hyman, J. Comp. Phys., Vol. 50, 1983. 

13 H.C. Yee, “A comparative Study of Flux Limiters for Two-Dimensional Time-Accurate Calculations,” 
in preparation. 

14 P.L. Roe, J. Comp. Phys. Vol. 43, 1981. 

15 H.C. Yee and P. Kutler, NASA TM-85845, August, 1983. 

16 M. Mandella and D. Bershader, “Generation and Aerodynamic Interaction of Compressible Vortices,” 
in preparation. 

17 Y. J. Moon and H.C. Yee, “Numerical Simulation by TVD Schemes of Complex Shock Reflections from 
Airfoils at High Angle of Attack,” in preparation. 




2. Government Accession No. 


3. Recipient’s Catalog No, 


1. Report No. 

NASA TM-88325 

2. Government Accession No. 

3. Recipient’s Catalog No, 

4. Title and Subtitle 

NUMERICAL EXPERIMENTS WITH A SYMMETRIC HIGH- 
RESOLUTION SHOCK-CAPTURING SCHEME 

5. Report Date 

June 1986 

6. Performing Organiiation Code 

7. Author|f) 

H. C. Yee 

8. Performing Organization Report No. 

86316 

10. Work Unit No. 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 
505-31-01-01-00-21 

9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 

15. Supplementary Notes 

Point of Contact: H. C. Yee, Ames Research Center, M/S 202 A- 1 

Moffett Field, CA 94035 (415) 694-5548 or FTS 464-5548 

16. Abstract 

Characteristic-based explicit and implicit total variation diminishing 
(TVD) schemes for the two-dimensional compressible Euler equations have 
recently been developed. This is a generalization of recent work of Roe 


and Davis to a wider class of symmetric (non-upwind) TVD schemes other than 
Lax-Wendroff . Roe and Davis's schemes can be viewed as a subset of the 
class of explicit methods. The main properties of the present class of 
schemes are that they can be implicit, and, when steady-state calculations 
are sought, the numerical solution is independent of the time step. In a 
recent paper, a comparison of a linearized form of the present implicit 
symmetric TVD scheme with an implicit upwind TVD scheme originally 
developed by Harten and modified by Yee was given. The results favored the 
symmetric method. It was found that the symmetric method is just as accu- 
rate as the upwind method while requiring less computational effort. Cur- 
rently, more numerical experiments are being conducted on time-accurate 
calculations and on the effect of grid topology, numerical boundary condi- 
tion procedures, and different flow conditions on the behavior of the 
method for steady-state applications. The purpose of this paper is to 
report our experiences with this type of scheme and give some basic guide- 
lines for using the scheme. 


17. Key Worth (Suggwted by AulhorUI) Numerical method, IB. Diltribution Statement 

Finite difference method, Computational Unlimited 

fluid dynamics, System of hyperbolic 
conservation laws, Weak solutions, 

Shock capturing, Conservative differ- 



Sub: 

ect Cateeorv 

- 64 

19. Security Oassif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price" 

Unclassified 

Unclassified 

10 

A02 


'For fate by the National Technical Information Service, Springfield. Virginia 22161 















